home *** CD-ROM | disk | FTP | other *** search
/ WINMX Assorted Textfiles / Ebooks.tar / Text - Mathematics - Numerical Mathematics and Computing (F).zip / amrkad.f < prev    next >
Text File  |  2002-06-11  |  5KB  |  196 lines
  1. C
  2. C PAGE 408-409: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
  3. C
  4. C FILE: AMRKAD.FOR
  5. C
  6. C ADAPTIVE SCHEME FOR ADAMS-MOULTON METHOD FOR SYSTEMS OF ODE'S
  7. C  (AMRKAD,XPSYS,AMRK,AMSYS,RKSYS)
  8. C
  9.       DIMENSION XA(10)      
  10.       TA = 3.0 
  11.       TB = .5    
  12.       XA(1) =3. 
  13.       XA(2) =6. 
  14.       XA(3) = 3.
  15.       H = -.2     
  16.       ITMAX = 100 
  17.       EMIN = 1.0E - 8  
  18.       EMAX = 1.0E - 2      
  19.       HMIN = 1.0E - 4   
  20.       HMAX = 1.0 
  21.       PRINT 13,TA,H 
  22.       PRINT 15,(XA(I),I=1,3)
  23.       CALL AMRKAD(TA,XA,H,TB,ITMAX,EMIN,EMAX,HMIN,HMAX,IFLAG)       
  24.       PRINT 5,IFLAG 
  25.       STOP
  26.     5 FORMAT('0  IFLAG =',I2)       
  27.    13 FORMAT(2X,2F10.5)     
  28.    15 FORMAT(5X,7E15.8)     
  29.       END 
  30.   
  31.       SUBROUTINE XPSYS(N,X,F) 
  32.       DIMENSION X(10),F(10) 
  33.       N = 2       
  34.       X1 = X(1)   
  35.       X2 = X(2)   
  36.       F(1) = 1.0  
  37.       F(2) = 3.0*X2/X1 + 4.5*X1 - 13. 
  38.       RETURN      
  39.       END 
  40.  
  41.       SUBROUTINE AMRKAD(T,X,H,TB,ITMAX,EMIN,EMAX,HMIN,HMAX,IFLAG)   
  42.       DIMENSION X(10),Z(10,5),F(10,5) 
  43.       DATA  EPSI/1.0E-6/   
  44.       M = 1       
  45.       IRK = 3     
  46.       IAM = 0     
  47.       IFLAG = 3   
  48.       NSTEP = 0   
  49.       IF(ABS(TB-T) .LT. ABS(4.0*H))  H = 0.25*H 
  50.       CALL XPSYS(N,X,F(1,M))
  51.       PRINT 10,T,H
  52.       PRINT 11,(X(I),I=1,N) 
  53.       DO 1 I=1,N  
  54.    1  Z(I,M) = X(I) 
  55. C
  56.    2  DT = ABS(TB -T)       
  57.       IF(DT .GE. ABS(H))  GO TO 3   
  58.       IFLAG = 0   
  59.       IF(DT .LE. EPSI*AMAX1(ABS(TB),ABS(T)))  GO TO 8     
  60.       H = SIGN(DT,H)
  61.       IRK = 1     
  62.       IAM = 0     
  63.       IF(HMIN .LE. ABS(H) .AND. ABS(H) .LE. HMAX)  GO TO 3
  64.       IFLAG = 2   
  65.       GO TO 8     
  66. C
  67.    3  IF(IAM .EQ. 1)  GO TO 5 
  68.       DO 4 K=1,IRK
  69.       CALL RKSYS(N,M,T,Z,H,F)
  70.       NSTEP = NSTEP + 1     
  71.       PRINT 10,T,H
  72.    4  PRINT 11,(Z(I,M),I=1,N) 
  73.       IF(NSTEP .GT. ITMAX .OR. IRK .EQ. 1)  GO TO 8       
  74.    5  CALL AMSYS(N,M,T,Z,H,F,EST)     
  75.       PRINT 10,T,H,EST      
  76.       PRINT 11,(Z(I,M),I=1,N) 
  77.       NSTEP = NSTEP + 1     
  78. C       
  79.       IF(NSTEP .GT. ITMAX .OR. IFLAG .EQ. 0)  GO TO 8     
  80.       IAM = 1     
  81.       HSAVE = H   
  82.       IF(EST .LT. EMIN)  GO TO 6      
  83.       IF(EST .LE. EMAX)  GO TO 2      
  84.       H = 0.5*H   
  85.       GO TO 7     
  86.    6  IF(DT .LT. ABS(4.0*H))  GO TO 2 
  87.       H = 2.0*H   
  88.    7  M = 1 + MOD(M,5)      
  89.       T = T - 4.0*HSAVE     
  90.       IAM = 0     
  91.       IF(HMIN .LE. ABS(H) .AND. ABS(H) .LE. HMAX)  GO TO 2
  92.       IFLAG = 1   
  93. C       
  94.    8  DO 9 I=1,N  
  95.    9  X(I) = Z(I,M) 
  96.       RETURN      
  97.   10  FORMAT(2X,'T,H,EST:',3(E10.3,1X)) 
  98.   11  FORMAT(8X,'X:',5(5X,E22.14))    
  99.       END 
  100.   
  101.       SUBROUTINE AMRK(N,T,X,H,NSTEP)  
  102.       DIMENSION  X(10),F(10,5),Z(10,5)
  103.       M = 1       
  104.       PRINT 6,T,H 
  105.       PRINT 7,(X(I),I=1,N)  
  106.       DO 2 I=1,N  
  107.         Z(I,M) = X(I)       
  108.    2  CONTINUE    
  109.       DO 3 K = 1,3
  110.         CALL RKSYS(N,M,T,Z,H,F)       
  111.         PRINT 6,T,H 
  112.         PRINT 7,(Z(I,M),I = 1,N)      
  113.    3  CONTINUE
  114.       DO 4 K = 4,NSTEP      
  115.         CALL AMSYS(N,M,T,Z,H,F,EST)   
  116.         PRINT 6,T,H,EST     
  117.         PRINT 7,(Z(I,M),I = 1,N)      
  118.    4  CONTINUE
  119.       DO 5 I=1,N  
  120.         X(I) = Z(I,M)       
  121.    5  CONTINUE
  122.       RETURN      
  123.    6  FORMAT(2X,'T,H,EST:',3(1X,E10.3)) 
  124.    7  FORMAT(8X,'X:',5(5X,E22.14))
  125.       END 
  126.   
  127.       SUBROUTINE RKSYS(N,M,T,X,H,F)   
  128.       DIMENSION  X(10,5),XP(10),F(10,5),F2(10),F3(10),F4(10)
  129.       MP1 = 1 + MOD(M,5)    
  130.       H2 = 0.5*H  
  131.       CALL XPSYS(X(1,M),F(1,M))       
  132.       DO 2 I = 1,N
  133.         XP(I) = X(I,M) + H2*F(I,M)    
  134.    2  CONTINUE
  135.       CALL XPSYS(XP,F2)     
  136.       DO 3 I = 1,N
  137.         XP(I) = X(I,M) + H2*F2(I)     
  138.    3  CONTINUE
  139.       CALL XPSYS(XP,F3)     
  140.       DO 4 I = 1,N
  141.         XP(I) = X(I,M) + H*F3(I)      
  142.    4  CONTINUE
  143.       CALL XPSYS(XP,F4)     
  144.       DO 5 I = 1,N
  145.         X(I,MP1) = X(I,M) + H*(F(I,M) + 2.0*(F2(I) + F3(I)) + F4(I))/6.0      
  146.   5   CONTINUE
  147.       T = T + H   
  148.       M = MP1     
  149.       RETURN      
  150.       END 
  151.   
  152.       SUBROUTINE AMSYS(N,M,T,X,H,F,EST) 
  153.       DIMENSION  X(10,5),XP(10),F(10,5),SUM(10),CP(4),CC(4) 
  154.       DATA  CP/55.0,-59.0,37.0,-9.0/  
  155.       DATA  CC/ 9.0, 19.0,-5.0, 1.0/  
  156.       MP1 = 1 + MOD(M,5)    
  157.       CALL XPSYS(X(1,M),F(1,M))       
  158.       DO 2 I = 1,N
  159.         SUM(I) = 0.0
  160.    2  CONTINUE
  161.       DO 4 K = 1,4
  162.         J = 1 + MOD(M-K+5,5)
  163.         DO 3 I = 1,N
  164.           SUM(I) = SUM(I) + CP(K)*F(I,J)
  165.    3    CONTINUE  
  166.    4  CONTINUE
  167.       DO 5 I = 1,N
  168.         XP(I) = X(I,M) + H*SUM(I)/24.0
  169.    5  CONTINUE    
  170.       CALL XPSYS(XP,F(1,MP1)) 
  171.       DO 6 I = 1,N
  172.         SUM(I) = 0.0
  173.    6  CONTINUE    
  174.       DO 8 K = 1,4
  175.         J = 1 + MOD(MP1-K+5,5)
  176.         DO 7 I = 1,N
  177.           SUM(I) = SUM(I) + CC(K)*F(I,J)
  178.    7    CONTINUE  
  179.    8  CONTINUE
  180.       DO 9 I = 1,N
  181.         X(I,MP1) = X(I,M) + H*SUM(I)/24.0       
  182.    9  CONTINUE    
  183.       T = T + H   
  184.       M = MP1     
  185.       DLTMAX = 0.0
  186.       DO 10 I = 1,N 
  187.         DLT = ABS( X(I,M) - XP(I) )   
  188.         IF(DLT .GT. DLTMAX) THEN      
  189.           DLTMAX = DLT      
  190.           IDLT = I
  191.         END IF    
  192.   10  CONTINUE    
  193.       EST = 19.0*DLTMAX/( 270.0*ABS(X(IDLT,M)) )
  194.       RETURN      
  195.       END 
  196.